Introduction

Sustainability of rice (Oryza sativa L.) production is highly dependent on the available water resources. Increasingly limited resources of fresh water are likely to become a major factor limiting the optimal growth of plants (Tietenberg and Lewis 2009). Many plant-breeding programs focus on producing new crop varieties with traits that will help enhance food security. Most modern rice cultivars are susceptible to abiotic stresses (Swamy et al. 2017). Rice uses a lot of water for growth; on average, 2500 liter of water is needed to produce 1 kg of rice (Bouman 2009). However, fresh water is becoming limited in various regions of the world due to global climate change (Malaysian Meteorological Department 2017). Thus, it is important to develop new rice lines that can tolerate limited water availability while still producing high yields.

The year 2016 was the hottest year on record in Malaysia, with an average temperature of 27.66°C (Malaysian Meteorological Department 2017). This year was even hotter than 1998, when the most severe drought on record struck Malaysia. In 2016, drought caused by a super El-Niño event occurred throughout Malaysia between February and May (Malaysian Meteorological Department, 2017). The 2016 El-Niño event created considerable anxiety among rice farmers, who needed to adjust their planting schedules to avoid the worst effects of the drought. Rice farmers in granary areas are less affected by drought because they can acquire an adequate water supply, but most farmers outside the main granary areas where access to water is limited and reliant solely on rain for irrigation.

The development of drought-tolerant rice cultivars is one approach for helping farmers reduce the impact of drought and obtain better yields. However, until 2018, no drought-tolerant cultivars have been released in Malaysia despite numerous cultivars being released in neighboring countries. In 2013, the rice cultivar MR1A was released for cultivation. This cultivar can be grown under aerobic conditions and requires less water than typical varieties, but its yield is low; 2.0 to 3.0 tons per hectare (Othman et al. 2014). Collaborative research by Universiti Kebangsaan Malaysia (UKM) and the International Rice Research Institute (IRRI) developed drought-tolerant lines by pyramiding drought yield QTLs (qDTY) into local cultivars MR219 and MRQ74. These two drought-tolerant rice lines showed better yields than local cultivars under drought stress trials (Shamsudin et al. 2016a, b). It is hoped that the pyramided lines (PLs) can be used by farmers to mitigate the disastrous effects of droughts.

However, yield performance of a genotype can vary from one environment to another due to GEI (Fasahat et al. 2014). Genes controlling yield might be expressed differently in a different environment, making one genotype superior to others in any given environment. A genotype is stable when its yield is consistent among environments. Therefore, evaluating genotypic stability is important prior to releasing a cultivar to farmers so they can be assured of obtaining similar yields in different locations.

According to Huehn (1990), parametric procedure is a good attribute when it relies on certain statistical assumption such as interaction effects and normal distribution errors. A drawback with this approach is that the statistical assumption is sensitive to the significance of variances and variance-related measures. Here is where non-parametric measures provide an alternative approach since the procedure is conducted without the underlying specific assumptions. The non-parametric methods are based on the ranks of genotypes in different environments. The genotypes are considered stable if their rank is similar despite the environmental differences. Regardless of merits and demerits from both approaches, each of them can complement and supplement each other, hence creating a better picture of interaction for GEI interpretation (Dehghani et al. 2016).

Researchers frequently use genotype and genotype-by-environment interaction (GGE) models and additive main effects and multiplicative interaction (AMMI) models to obtain yield trial data on genetic crosses (Gauch 2006). The concept of using GGE biplots, as proposed by Yan et al. (2000) is that genotype (G) and genotype-by-environment interaction (GE) must be considered simultaneously when making selection decisions. The biplot data can be used to determine crossover interactions or rank changes in yields of cultivars under various environmental conditions.

Choices for GGE biplot scaling include genotype-focused scaling, environment-focused scaling, symmetric scaling, and equal-space scaling, all of which display a “which-won-where” pattern. Although none of these scaling approaches is perfect, each complements the others. It is believed that the use of GGE biplots can identify highly adapted and phenotypically stable lines across a wide range of environments (Oladosu et al. 2017). In contrast, the AMMI model proposed (Gauch 1992) combines a univariate (ANOVA) to assess the main effect of the genotypes and environment together with a multivariate technique; usually principal component analysis (PCA) to assess genotype-by-environment interactions. The principal components of a PCA usually represent the response of genotypes that are proportional and not proportional to the environment. AMMI models are usually called AMMI(n), with n defined as the number of components used to study the interactions. The use of these two models has been critically reviewed and compared (Gauch 2006; Yan et al. 2007; Gauch et al. 2008). Numerous studies have used the GGE biplot method to study the genotypic stability of varieties (Akter et al. 2015; Balakrishnan et al. 2016; Shahriari et al. 2018). We believe that the use of the GGE biplot approach is more informative than the AMMI approach, despite the criticisms conveyed by Oladosu et al. (2017). The present study uses the GGE model, univariate and non-parametric approaches to analyze the performance and genotypic stability of drought-tolerant PLs under drought and water-limited environments across Malaysia.

Materials and Methods

Plant materials, location and management practices

Table 1: List of rice genotypes tested and the information

Designation	Type	Source
IR 98010-134-4-1-2-1-1	MRQ74 PL	UKM
IR 99784-156-137-1-3-1-1	MR219 PL	UKM
IR 99784-226-335-1-2-1-1	MR219 PL	UKM
IR 99784-226-335-1-5-1-1	MR219 PL	UKM
IR 99784-255-68-1-7-1-1	MR219 PL	UKM
IR 99784-255-7-2-5-1-1	MR219 PL	UKM
IR 99784-255-91-1-1-1-1	MR219 PL	UKM
MR219	Stable check (modern cultivar)	MARDI
MRQ74	Check (modern cultivar)	MARDI
IR 77298-14-1-2-10	Check (drought tolerant line)	IRRI

Table 2: Description and characteristics of the environments tested in this study

Environment	Location	Year	Irrigation regime	Range of soil moisture content (%)	Mean soil water potential (kPa)	Maximum temperature (°C)	Minimum temperature (°C)	Mean temperature (°C)	Mean monthly rainfall (mm)	Soil series
TC1	Teluk Chengai	2015	Normal	65-68	-4	30.2	25.0	28.1	260.60	Chengai
TC1RS	Teluk Chengai	2015	Drought	24-27	-60	30.2	25.0	28.1	260.60	Chengai
TC2	Teluk Chengai	2016	Normal	62-68	-8	28.7	24.5	27.3	142.45	Chengai
TC2RS	Teluk Chengai	2016	Water-limited	45-51	-30	28.7	24.5	27.3	142.45	Chengai
TC3	Teluk Chengai	2017	Normal	69-71	-3	29.8	24.5	27.2	248.15	Chengai
TC3RS	Teluk Chengai	2017	Water-limited	48-52	-35	29.8	24.5	27.2	248.15	Chengai
SR	Sawah Ring	2017	Normal	67-72	-10	31.0	23.0	27.0	150.00	Sedu
SRRS	Sawah Ring	2017	Drought	26-28	-65	31.0	23.0	27.0	150.00	Sedu
PB	Parit Buntar	2017	Normal	71-73	-4	28.2	24	26.7	293.73	Selangor
PBRS	Parit Buntar	2017	Drought	21-23	-75	28.2	24	26.7	293.73	Selangor
BM	Kampung Bukit Merah	2016	Normal	65-70	-7	29.2	24.9	27.2	181.30	Kranji
BMRS	Kampung Bukit Merah	2016	Drought	26-30	-79	29.2	24.9	27.2	181.30	Kranji
TS	Parit Buntar	2017	Normal	68-71	-8	25.0	29.9	27.9	117.20	Selangor

Six drought-tolerant pyramided lines (PLs) of MR219 and one of MRQ74 selected from a previous advanced yield trial (AYT) (Shamsudin et al. 2016a, b; Ikmal et al. 2018, 2019) together with three checks namely IR 77298-14-2-10 (drought tolerant line) and the two recipient parents; MR219 and MRQ74 (Table 1). These drought-tolerant PLs generated from crosses between donors of drought yield QTL (qDTY) from the International Rice Research Institute (IRRI), Philippines and the recipient parents, MR219 and MRQ74 from Malaysia. These PLs have different combinations of qDTYs viz., qDTY_2.2, qDTY_3.1, and qDTY_12.1. Table 2 shows the 13 environments used to conduct this study which was carried out in 2015, 2016 and 2017. Environments were defined as the combination of location, treatment and years. Randomized complete block design (RCBD) with three replications was used in this study. Plot size, 2 m × 5 m with planting spacing of 25 cm between rows and hills was used. Rice cultivation guidelines given by the Department of Agriculture Malaysia (DOA) were used for field cultural practices. Fertilizer (17 N: 20 P: 10 K) was given at two separate occasions; firstly, during the early growth (seven days after germination) at the rate of 140 kg per ha and secondly after 50 days of germination at the rate of 100 kg per ha. Urea fertilizer was applied during the active tillering period at the rate of 80 kg per ha. Chemical fungicides and insecticides were sprayed to the field to control diseases and pests. Standing water in the plot was only allowed for 30 days after transplanting and drained on the 31st day to create drought stress or water-limited condition. If rain happened to occur, the drains at the plot were opened to allow water to flow out from the plot to ensure the low soil moisture. In TS (Parit Buntar), the drought evaluation was carried out in a concrete containing soil and equipped with out-flow channels made from polyvinyl chloride (PVC) pipes to drain out excessive water to create drought condition. For water-limited condition, water was supplied at the field capacity instead of 5 cm standing water above the soil surface for normal condition.

Phenotyping

Grain yield was obtained after grains were harvested, dried, weighed and adjusted to 14% moisture content from ten inner rows, leaving 50 cm at both ends per plot (or two hills) as border area of the experimental plot. The final weights were converted to kilogram per hectare (kg ha^-1). The Standard Evaluation System (SES) for Rice (IRRI 2013) used as a guide for all other yield related traits such as the number of panicle (NP), spikelet per panicle (SPP), filled spikelets (FS) and thousand-grain weight (TGW). NP, SPP, FS and TGW were recorded from trials in Parit Buntar, Teluk Chengai and Bukit Merah.

Data analysis

Combined analysis of variance (ANOVA) was computed for genotype (G), environment (E) and GE interaction for grain yield using Statistical Tool for Agricultural Research (STAR) version 2.0.1. Rstudio used to produce GGE biplot using the GGEBiplotGUI package (Frutos et al. 2014). Stability measures used in this study viz., linear regression coefficient (b_i), Shukla stability variances (σ_i²) (Shukla 1972), deviation from regression (S²_d) (Eberhart and Russell 1966) and Wricke’s ecovalence (W_i²) (Wricke 1962). Kang’s stability statistics (YS_i) (Magari and Kang 1993) were computed in R studio using Agricolae package. Non-parametric measures of stability (Nassar and Huhn 1987; Kang 1988; Thennarasu 1995) also computed in R studio using Phenability package. The relationship between all stability measures calculated in R Studio using Spearman rank correlation and visualized using Corrplot package.

Results

All locations selected as the test environments received medium amount of monthly rainfall that was between 100 – 300 mm. The period of June to July is classified as the driest month in most states of Malaysia while November, December and January are months with the maximum amount of rainfall by Malaysian Meteorological Department. Parit Buntar recorded the lowest and the highest amount of mean monthly rainfall during the planting seasons (Table 2). Drought condition is achieved when the soil moisture content dropped below 30% and at the permanent wilting point while water-limited condition is when the soil moisture content between 50–60%.

Combined analysis of variance

Table 3: Combined analysis of variance for grain yield in 13 environments

Source	DF	Sum of Square	Mean Square	F-Value	Pr (> F)	%SS
Environment (E)	12	750020796.13	62501733.01	23.30	***	72.00
Genotype (G)	9	75040580.40	8337842.27	7.74	***	7.20
Genotype × Environment (GE)	108	216615870.36	2005702.50	1.86	***	20.79
Replicate (Environment)	26	69746412.48	2682554.33	2.49	***
Pooled Error	234	252207608.78	1077810.29
Total	389	1363631268.15

Percentage sum of square (%SS)

Table 4: Mean values of grain yield, number of panicle, length of panicle, spikelet per panicle, filled spikelets and thousand-grain weight for each genotype combined across environments

Genotype	Grain yield (kg ha^-1)	Number of panicle	Length of panicle (cm)	Spikelet per panicle	Filled spikelets	Thousand-grain weight (g)
IR 98010-134-4-1-2-1-1	5393.99bcd	17.97 a	25.72a	184.25a	164.50a	23.94a
IR 99784-156-137-1-3-1-1	6328.79a	14.26 c	25.44a	179.67a	161.33a	27.30a
IR 99784-226-335-1-2-1-1	6024.26abc	15.13 bc	25.44a	169.25a	158.00a	27.20a
IR 99784-226-335-1-5-1-1	6388.31a	16.38 ab	24.39a	153.67a	142.33a	28.23a
IR 99784-255-68-1-7-1-1	6086.41ab	15.90 bc	24.14a	180.67a	162.50a	26.63a
IR 99784-255-7-2-5-1-1	6256.98a	15.15 bc	25.23a	165.92a	154.17a	27.24a
IR 99784-255-91-1-1-1-1	6200.37a	16.49 ab	25.23a	168.83a	159.58a	26.12a
MR219	5705.24abcd	15.41 bc	25.33a	173.08a	160.33a	26.87a
MRQ74	5287.13cd	16.95 ab	25.32a	162.67a	148.83a	25.70a
IR 77298-14-1-2-10	5146.55d	16.23 ab	25.48a	162.83a	148.17a	27.05a
Mean	5881.80	15.99	25.17	170.13	155.97	26.63
SE	94.81	0.21	0.19	3.24	3.00	0.30
CV (%)	17.65	26.34	8.41	20.87	21.07	12.33
Min	1293.22	7.00	20.40	73.00	67.00	18.80
Max	11720.00	34.00	30.60	289.00	287.00	35.00

Coefficient of variation (CV), standard error (SE), minimum (min), maximum (max). Mean values with different letter are significantly different (Tukey’s HSD, P < 0.05)

Combined analysis of variance revealed that genotype (G), environment (E) and genotype-by-environment (GE) interaction for grain yield were highly significant (P < 0.001) (Table 3). We found that G explained 7.20% of the total variance, E explained 72.00% of the total variance, and GE explained 20.79% of the total variance.

The mean values for genotype comparison

Table 4 shows the grand mean values of grain yield and the yield related traits for each genotype across environments. IR 99784-226-335-1-5-1-1 has the highest grain yield (6388.31 kg ha^-1) followed by IR 99784-156-137-1-3-1-1 (6328.79 kg ha^-1) while IR 77298-14-1-2-10 has the lowest grain yield (5146.55 kg ha^-1). The best genotype showed 21 and 24% higher yield than MR219 and MRQ74 respectively. The drought-tolerant MR219 PLs produced 319.02–683.07 kg ha⁻¹ higher yield than MR219 (G8), whereas MRQ74 PL produced 106.86 kg ha⁻¹ higher yield than MRQ74 (G9). Mean number of panicle (NP) for all genotypes were significantly different from each other across environment with IR 98010-134-4-1-2-1-1 recorded the highest while IR 99784-156-137-1-3-1-1 recorded the lowest NP. Length of panicle (LP) for all genotypes was more than 25.00 cm except IR 99784-226-335-1-5-1-1 and IR 99784-255-68-1-7-1-1. Meanwhile, number of spikelet per panicle (SPP) for every genotypes were more than 150 with IR 98010-134-4-1-2-1-1 had the highest SPP. The same result also can be seen for number of filled spikelet per panicle (FS), where genotypes with high number of SPP will have high number of FS. IR 99784-226-335-1-5-1-1 recorded the highest thousand-grain weight (TGW) while IR 98010-134-4-1-2-1-1 had the lowest TWG. The mean values of grain yield and the other yield contributing traits for each environment is provided in Supplementary Tables S1, S2 and S3.

GGE biplot analysis

The GGE biplot explained 66.09% (PC1=47.13%, PC2=18.96%) of the total variation that was related to G and GE interaction. The what-won-where pattern for grain yield was shown (Fig. 1A). The what-won-where polygon pattern of GGE biplot was constructed in such a way that all the tested genotypes were contained within the polygon. The straight line originating from the center of the biplot divides the polygon into different sections. The genotypes at the vertex or edge of the polygon in a sector represent the winning genotype for entire environments in that sector. Also, the lines that originate from the biplot center perpendicular to the polygon represent a hypothetical environment in which the two genotypes that represent the two sides of the polygon are said to perform equally. IR 99784-156-137-1-3-1-1 is the best genotype in the environments TC1, TC1RS, PB, PBRS, SR, SRRS, TC2, and TC3RS. Meanwhile, IR 99784-226-335-1-5-1-1 is the winning genotype in TC3, TS, and TC2RS. IR 98010-134-4-1-2-1-1 is the best in BM and BMRS. IR 99784-156-137-1-3-1-1 and IR 99784-226-335-1-5-1-1 performed the same in environments TC2 and TC3RS while IR 99784-226-335-1-5-1-1 and IR 98010-134-4-1-2-1-1 performed similarly in environment TC2RS and TS.

Fig. 1: (A) which won where biplot for grain yield showing different genotypes winning in different environments during the 2016 to 2017 trials. (B) Mean versus stability biplot for grain yield showing the stability and performance of each genotypes. (C) Ranking of genotypes compared to the “ideal genotypes” located in the innermost concentric circle as pointed by the arrowhead. (D) Discriminativeness versus representativeness biplot showing the discriminating ability and the representativeness of environments as measured by the angles to the AEC abscissa. G1: IR 98010-134-4-1-2-1-1, G2: IR 99784-156-137-1-3-1-1, G3: IR 99784-226-335-1-2-1-1, G4: IR 99784-226-335-1-5-1-1, G5: IR 99784-255-68-1-7-1-1, G6: IR 99784-255-7-2-5-1-1, G7: IR 99784-255-91-1-1-1-1, G8: MR219, G9: MRQ74, G10: IR 77298-14-1-2-10

Fig. 2: (A) Ranking of environments compared to the “ideal environment” located in the innermost concentric circle as pointed by the arrowhead. (B) Relationship among environments biplot showing the angles between vectors of environments to each other

Table 5: Stability statistics parameters

Genotype	Mean	bi	S^2d	σi^2	Wi^2	YSi	S^1	S^2	S^3	S^6	N^1	N^2	N^3	N^4	KRS
IR 98010-134-4-1-2-1-1	5393.99	0.64	3651191.3**	4306655.5**	43750736	-7	0.04	10.69	25.44	7.41	2.77	0.35	0.46	0.01	10
IR 99784-156-137-1-3-1-1	6328.79	1.26	2059493.4*	2386139.5*	25313782	7+	0.05	11.03	16.40	4.96	2.77	0.69	0.72	0.01	12
IR 99784-226-335-1-2-1-1	6024.26	1.05	988267.0ns	882327.2ns	10877184	6+	0.08	6.91	11.58	4.47	2.15	0.43	0.49	0.01	11
IR 99784-226-335-1-5-1-1	6388.31	1.23	995947.1ns	1313325.1ns	15014764	12+	0.10	9.27	10.44	3.33	2.46	0.82	0.81	0.03	10
IR 99784-255-68-1-7-1-1	6086.41	1.10	1975828.8*	1844699.3ns	20115956	5+	0.09	11.58	17.70	5.26	2.62	0.65	0.69	0.02	12
IR 99784-255-7-2-5-1-1	6256.98	1.20	1793834.1ns	1903810.7ns	20683426	7+	0.06	6.56	7.93	3.36	2.15	0.72	0.60	0.02	11
IR 99784-255-91-1-1-1-1	6200.37	0.91	975675.4ns	915937.3ns	11199841	8+	0.04	4.50	10.64	3.80	1.69	0.34	0.43	0.01	7
MR219	5705.24	1.13	490829.8ns	547590.1ns	7663708	3+	0.03	4.76	14.83	5.23	1.46	0.21	0.31	0.00	13
MRQ74	5287.13	0.65	1725149.8ns	2517108.1**	26571080	-8	0.01	9.44	23.69	8.35	2.62	0.29	0.42	0.00	10
IR 77298-14-1-2-10	5146.55	0.82	3532576.0**	3439432.3**	35425393	-9	0.01	13.73	27.71	8.86	3.31	0.41	0.47	0.00	14

Non-significant (ns). *significant at P < 0.05, **significant at P < 0.01

Fig. 3: Graphical correlation matrix showing the relationship among tested environments and the significant levels. *** P < 0.001, ** P < 0.01, * P < 0.05

The mean versus stability of the tested genotypes shows that the green line with a single arrowhead is the average environment coordinate (AEC) abscissa (Fig. 1B). The green lines perpendicular to the AEC are the AEC ordinate in which a particular genotype is considered less stable when its projection line is longer. The vertical line passing through the origin and is perpendicular to the AEC abscissa divides the genotypes into higher than the overall mean performance and lower than the overall mean performance across environments. As shown that, IR 99784-156-137-1-3-1-1, IR 99784-226-335-1-2-1-1, IR 99784-226-335-1-5-1-1, IR 99784-255-68-1-7-1-1, IR 99784-255-7-2-5-1-1, and IR 99784-255-91-1-1-1-1 are the genotypes with higher mean grain yields than the overall mean grain yield while IR 98010-134-4-1-2-1-1, MR219, MRQ74, IR 77298-14-1-2-10 are the genotypes with lower mean grain yields than the overall mean grain yield (Fig. 1B). The arrowhead is pointing towards the genotypes with higher mean grain yields and consequently ranked the genotypes with respect to the grain yield. IR 99784-226-335-1-5-1-1 is the most stable and high yielding followed by IR 99784-226-335-1-5-1-1 and IR 99784-255-7-2-5-1-1. Although IR 99784-156-137-1-3-1-1 won in more environments than IR 99784-226-335-1-5-1-1, its stability is lower as indicated by the longer projection of IR 99784-156-137-1-3-1-1 than IR 99784-226-335-1-5-1-1. MR219 had almost similar stability to IR 99784-255-7-2-5-1-1 but grain yield of MR219 is much lower than IR 99784-255-7-2-5-1-1. IR 98010-134-4-1-2-1-1 is the least stable genotype evidenced by the longest projection from the AEC abscissa.

The concentric circles are used to illustrate the distance between the genotypes and the “ideal genotype” that is located at the center or the innermost circle (Fig. 1C). An “ideal genotype” should both be high yielding and have great stability. IR 99784-226-335-1-5-1-1 had the shortest distance from the innermost circle followed by IR 99784-255-7-2-5-1-1 and IR 99784-156-137-1-3-1-1. Three genotypes namely IR 98010-134-4-1-2-1-1, MRQ74, and IR 77298-14-1-2-10 had the furthest distance from the innermost circle. The ranking of the genotypes for grain yield is IR 99784-226-335-1-5-1-1> IR 99784-255-7-2-5-1-1> IR 99784-156-137-1-3-1-1> IR 99784-255-91-1-1-1-1> IR 99784-226-335-1-2-1-1> IR 99784-255-68-1-7-1-1> MR219> MRQ74> IR 98010-134-4-1-2-1-1> IR 77298-14-1-2-10.

Discriminativeness and the representativeness of the environments are illustrated in Fig. 1D. The former ability is measured by the length of the environmental markers' vector while the latter is measured by the size of the angle of the environmental marker's vector to the AEC abscissa. TC1 has the longest vector from the origin of the biplot followed by TC1RS and TS, while SRRS has the shortest vector. The angles between the vectors of TC3RS and TC2 are the top two smallest to the AEC abscissa compared to other environments. BMRS has the largest angle to the AEC abscissa. Fig. 2A shows the ranking of environments compared to the ‘ideal environment’ where TC1RS and TC3 are the closest environments to the ‘ideal environment’.

The relationship between test environments have been shown that fifty-two environment combinations were positively correlated, and the remaining 26 environment combinations were negatively correlated (Fig. 2B). The significantly positively correlated environment combinations are TS/TC2RS, PB/TC1, PB/TC3RS, PB/TC1RS, TC1/TC1RS, TC1/TC2, TC3/TC1RS, TC3/TC2, TC3RS/TC1RS, TC3RS/TC2, and TC1RS/TC2. The graphical correlation matrix between the test environments showed that positively correlated environments are closely related and any genotypes tested in that environments will produce almost similar grain yield but the opposite will happen if the environments are negatively correlated to each other (Fig. 3).

The linear regression coefficient (b_i) ranged from 0.64 to 1.26. IR 99784-226-335-1-2-1-1 recorded the best value of b_i, 1.05 that is the closest to 1.00 at P < 0.05, followed by IR 99784-255-91-1-1-1-1 with the b_i of 0.91 at P < 0.05. MR219, IR 99784-255-91-1-1-1-1, and IR 99784-226-335-1-2-1-1 are the top three genotypes with the lowest values of S²_d compared to the other genotypes (Table 5). Their values of deviation form regression (S²_d) were not significantly different from zero. IR 98010-134-4-1-2-1-1 showed the highest S²_d value was significantly different from zero. According to the values of Shukla stability variance (σ_i²) and Wricke’s ecovalence (W_i²), MR219 has the lowest values, followed by IR 99784-226-335-1-2-1-1 and IR 99784-255-91-1-1-1-1, while IR 98010-134-4-1-2-1-1 has the highest values. For the yield stability statistic (YS_i), seven genotypes were marked with the “+” signs, are IR 99784-156-137-1-3-1-1, IR 99784-226-335-1-2-1-1, IR 99784-226-335-1-5-1-1, IR 99784-255-68-1-7-1-1, IR 99784-255-7-2-5-1-1, IR 99784-255-91-1-1-1-1, and MR219. IR 99784-226-335-1-5-1-1 has the highest value of YS_i. As shown in Table 5, MRQ74, and IR 77298-14-1-2-10 have the lowest values of the non-parametric stability measure of the mean of the absolute rank differences of a genotype over the n environments (S¹) compared to the other genotypes. Meanwhile, IR 99784-255-91-1-1-1-1, MR219 and IR 99784-255-7-2-5-1-1 had the lowest values of the variance among the ranks over the k environments (S²). For the other two non-parametric statistics, the sum of the absolute deviations (S³) and the relative sum of squares of rank for each genotype (S⁶) showed that IR 99784-226-335-1-5-1-1, IR 99784-255-7-2-5-1-1, and IR 99784-255-91-1-1-1-1 were the top three genotypes with the lowest values, whereas IR 77298-14-1-2-10 recorded the highest values. The values of N¹ showed that IR 99784-226-335-1-2-1-1 and IR 99784-255-7-2-5-1-1 are the lowest ranked genotypes. MR219 has the lowest value of N² while IR 99784-255-68-1-7-1-1 has the highest value. For N³, MR219 has the lowest value while the highest value recorded by IR 99784-226-335-1-5-1-1. Three genotypes recorded 0.00 values for N⁴ despite the unclear differences of the N⁴ values for all genotypes. It was also found that IR 99784-255-91-1-1-1-1 is the genotype with the lowest value of the Kang’s rank sum (KRS) stability statistics, while IR 77298-14-1-2-10 has the highest value of 14.00.

Discussion

The ten genotypes, including G8 (MR219), which is a stable mega-variety (Zainuddin et al. 2012) and G9 (MRQ74), which is a high-quality, specialty rice cultivar (Asfaliza et al. 2008; Suhaimee et al. 2009) were evaluated in present study. In partitioning of total variation, the large percentage of variation (due to GE interaction than to G) suggests that there were large differences in the genotype’s performance across environments. The existence of a different mega-environment is inferred where there is a significant amount of variation in GE interactions; that is, for each mega-environment, a different genotype has the best yield (Yan and Kang 2002).

A previous study also reported a significant effect of the environment on mutant rice genotypes (Oladosu et al. 2017). Changes in the relative rankings of genotypes across environments suggest that even though yield depends on genetics and environmental factors may also play a significant role in modifying yields (Oladosu et al. 2017). Therefore, it is important to test genotypes from targeted environments to select the best-performing genotypes for any given environment or across all environments examined. A multilocation trial involves several rice genotypes for testing, but GE interactions increase the difficulty of selecting the best genotypes for release to farmers. Therefore, to benefit the most farmers, breeders must make an extra effort when selecting the best genotype for all environments.

The which-won-where pattern is based on relative genotypic stabilities and mean performances of genotypes in tested environments. This pattern has also the ability to explain the presence or absence of crossover GE interactions in explaining the potential existence of various mega-environments (Yan and Rajcan 2002). Yan et al. (2000) stated that crossover and non-crossover types of GE interactions constitute MET data, whereby the former term shows a shift in yield ranking, whereas the latter term indicates constant yield performance. If a single variety wins in all environments tested, then no crossover GE interaction would be detected, and a single mega-environment exists. The genotypes located within polygons are judged as being less responsive to environmental conditions than genotypes located at the corners or vertices (Yan et al. 2007). Genotypes in vertices perform poorly across all environments if no environmental marker falls within its sector. The test environments with different winning genotypes (vertex genotypes) are defined as the mega-environments (Sserumaga et al. 2015).

The qDTYs showed a positive effect on the enhancement yields under drought or water-limited conditions, as evidenced by the higher mean values of grain yield than MR219 and MRQ74 (for all PLs), even though the yield advantage of IR 98010-134-4-1-2-1-1 over MRQ74 was only 106.86 kg ha⁻¹. Previous studies also reported that lines with qDTYs performed better than parents without qDTYs (Shamsudin et al. 2016a, b; Ikmal et al. 2018, 2019). Interestingly, all genotypes (excluding IR 99784-226-335-1-5-1-1) recorded at least one better grain yielding genotype under water-limited conditions than under normal-irrigation conditions (Supplementary Table S1). Better grain yields under water-limited conditions have been previously reported (Zhang et al. 2009; Poli et al. 2018). One possible reason for this disparity might be the presence of qDTY_3.1 in genotypes that reportedly cause lower yields under irrigated conditions (Venuprasad et al. 2009; Dixit et al. 2014, 2017; Ikmal et al. 2018). Yield depends on the severity of water stress and the inherent genotypic differences where limited irrigation improves root growth, photosynthetic rate, and root oxidation activities (Poli et al. 2018).

IR 99784-226-335-1-5-1-1, IR 99784-255-7-2-5-1-1, and IR 99784-255-91-1-1-1-1 showed better genotypic stability, evidenced by their shorter projections than MR219, which is known to be stable and has been grown successfully in 90% of Malaysia’s rice-growing regions for more than 20 seasons (Zainuddin et al. 2012). Grain yield is associated with many yield-related traits, such as TGW, NP, LP, SPP and FS. The genotype with the highest grain yield (IR 99784-226-335-1-5-1-1) had a high NP and thousand-grain weight; even though its LP, SPP, and FS were lower than the other genotypes we tested (Supplementary Table S3). Furthermore, the seed width of genotype IR 99784-226-335-1-5-1-1 is broader than the seeds of the other genotypes (Supplementary Fig. S1 and S2), which contribute to a higher TGW. Previous study also reported that mutant rice with high tiller numbers produces higher yields in stressed environments (Poli et al. 2018). Because the number of panicles (or tillers) is one determinant of yield, we suggest that this trait be used in selecting stable and high-yielding rice genotypes in Malaysia (Table S2).

The length of an environmental vector from its biplot origin is used to determine the discriminating ability of a given environment. Longer vectors represent environments that have a higher ability to discriminate among genotypes. This means that genotypes associated with short vectors will perform similarly in those environments. The angle between vectors and the AEC abscissa for environments is used to determine the representativeness of those environments to the mega-environment. An environment is more like the mega-environment than to other environments when its angle to the AEC abscissa is smaller than the angles of other environments to the abscissa. Based on vector lengths and angles, test environments can be grouped into three categories: (1) environments that give no or little information on differences among genotypes, (2) environments that are useful for selecting superior genotypes, and (3) environments that are useful for discarding low-performing or the least-stable genotypes. If a tested environment is highly representative of the mega-environment, the trial for that environment is sufficient for testing traits; thus, the cost of cultivars evaluation is reduced.

TC1 had the longest vector for grain yield, which means that it is the best environment to discriminate genotypes for selecting the best genotype (Fig. 1D). TC1 also had the smallest angle, which makes it the most similar environment to the mega-environment. Most drought environments had a shorter projection than typical irrigated environments. This situation shows that most of the drought-tolerant PLs we tested perform similarly, despite carrying different qDTYs. However, mild drought stress might be another reason that the yields of these genotypes were so similar. Severe stress conditions are needed to clearly identify the most drought-tolerant genotypes.

An environment can be identified as an “ideal environment” if it has a high ability to discriminate genotypes and representativeness. The ideal environment in present study was denoted by the arrowhead located at the innermost circle of biplot (Fig. 2A). If an ideal environment does not exist, the ideal environment in the biplot serves to indicate the similarity of a test environment to the ideal environment (Yan and Kang 2002; Oladosu et al. 2017). The environments located closest to the ideal environment in our biplot were the granary areas of Malaysia, which are considered the highest yielding and most productive rice cultivated areas in Malaysia.

The genotypes with b_i values that approached 1 (unity value) and S²_d of zero are considered the most stable. If a genotype has b_i < 1, it can adapt better to the unfavorable environmental conditions, while b_i > 1 shows that the genotype can be adapted to favorable or high yielding environment. IR 99784-226-335-1-2-1-1 and IR 99784-255-91-1-1-1-1 have the best b_i values that approached 1.00 and accompanied by high GY are considered the best genotype with greater stability based on this measure. On the other hand, the genotype with a low value of S²_d which is not significantly different from zero is considered more stable. Therefore, IR 99784-226-335-1-2-1-1, IR 99784-255-91-1-1-1-1, and MR219 are the most stable genotypes that met the required criteria. According to Shukla (1972) and Wricke (1962), the genotypes that have the lowest stability variance are considered as the most stable. Eberhart and Russell (1966) also stated that the desired variety is the one with high yield mean, b_i = 1.00 and as low as possible S²_d. According to the yield stability statistic (YS_i), the genotypes that were marked with the “+” sign have the YS_i values which are above the YS_i mean besides having greater stability across the environments. The other stability statistics which was developed by Thennarasu (1995), and Nassar and Huhn (1987) ranked the genotypes for stability from the lowest value to the highest value. The stability parameters developed by Thennarasu (1995) were based on the adjusted ranks of genotypes within each test environment. Genotypes with low S¹ values have low mean of the differences in rank over the 13 environments while genotypes with low S² have low variance among the rank over all tested environments. Meanwhile, genotypes with low S³ have low sum of all the absolute deviations of yield in all the environments. A low S⁶ value indicated that the genotype have low relative sum of squares of rank which means that the GY in every environment are almost the same. The KRS proposed by Kang (1988) assigned the lowest rank (1) to the most stable genotype, which has the highest yield and the lower stability variance (Vaezi et al. 2018).

Conclusion

In this study, environmental variation was the largest contributor to the variation in grain yield. Therefore, breeders need information on cultivar yield (specific to rice cultivating areas) to provide appropriate guidance to farmers on what cultivar to plant, particularly where yield problems occur and persist for long periods. It was found that the genotypes evaluated in this study grew best in specific environments, and so each environment had a different “winning” genotype. These winning genotypes should be recommended to farmers for cultivation in an environment-specific manner. Based on our results of genotypes with acceptable yields and stability, it is recommended that the following genotypes be included in further trials for cultivar release: IR 99784-156-137-1-3-1-1, IR 99784-226-335-1-5-1-1, IR 99784-255-7-2-5-1-1, and IR 99784-255-91-1-1-1-1.

Acknowledgements

We would like to thank the Department of Agriculture Malaysia (Paddy, Industrial Crops and Floriculture Division) for their field assistance and support in the research station. We are grateful to Universiti Kebangsaan Malaysia for research grants GUP-2017-033 and DCP-2017-004/2. We would also like to extend our gratitude to Dr. Arvind Kumar and Dr. B.P. Mallikarjuna Swamy from the International Rice Research Institute for their suggestions and recommendations. Thanks to Dr. Yusuff Oladosu from Institute of Tropical Agriculture, Universiti Putra Malaysia for helping in results interpretation and Prof. Manjit S. Kang for his critical comments and suggestions for improvement of the manuscript.

Author Contributions

AMI, AASN conducted the field experimental design and phenotyping, data analysis and manuscript preparation. TATNAR, ZPDE conducted phenotyping and analysis in two of the trial sites. ZSA helped in preparing the manuscript and RW revised the manuscript. All authors read and finalized the manuscript.

References

Akter A, MJ Hasan, MU Kulsum, MH Rahman, AK Paul, LF Lipi, S Akter (2015). Genotype × Environment interaction and yield stability analysis in hybrid rice (Oryza Sativa L.) by AMMI biplot. Bang Rice J 19:83–90

Asfaliza R, O Omar, S Abdullah, AR Bakar, M Harun (2008). MRQ74 – a new aromatic rice variety. Bull Teknol Tanaman 5:15–20

Balakrishnan D, D Subrahmanyam, J Badri, AK Raju, YV Rao, K Beerelli, S Mesapogu, M Surapaneni, R Ponnuswamy, G Padmavathi, VR Babu, S Neelamraju (2016). Genotype × Environment interactions of yield traits in backcross introgression lines derived from Oryza sativa cv. Swarna/Oryza nivara. Front Plant Sci 7; Article 1530

Bouman S (2009). How Much Water Does Rice Use?. Rice Today 8:28–29

Dehghani MR, MM Majidi, A Mirlohi, G Saeidi (2016). Study of genotype by environment interaction in tall fescue genotypes and their polycross progenies in Iran based on AMMI model analysis. Crop Past Sci 67:792–799

Dixit S, A Singh, N Sandhu, A Bhandari, P Vikram, A Kumar (2017). Combining drought and submergence tolerance in rice: Marker-assisted breeding and QTL combination effects. Mol Breed 37:143–154

Dixit S, A Singh, MT Sta Cruz, PT Maturan, M Amante, A Kumar (2014). Multiple major QTLs lead to stable yield performance of rice cultivars across variable drought intensities. BMC Genet 15; Article 16

Eberhart SA, WA Russell (1966). Stability parameters for comparing varieties. Crop Sci 6:36-40

Fasahat P, K Muhammad, A Abdullah, MAR Bhuiyan, MS Ngu, HG Gauch, R Wickneswari (2014). Genotype × Environment assessment for grain quality traits in rice. Commun Biometr Crop Sci 9:71–82

Frutos E, MP Galindo, V Leiva (2014). An interactive biplot implementation in R for modeling genotype-by-environment interaction. Stochast Environ Res Risk Assess 28:1629–1641

Gauch HG (2006). Statistical analysis of yield trials by AMMI and GGE. Crop Sci 46:1488–1500

Gauch HG (1992). Statistical Analysis of Regional Yield Trials: AMMI Analysis of Factorial Designs. Elsevier, Amsterdam, The Netherlands

Gauch HG, HP Piepho, P Annicchiarico (2008). Statistical analysis of yield trials by AMMI and GGE: Further considerations. Crop Sci 48:866–89-113

Huehn M (1990). Nonparameteric measures of phenotypic stability. part i. theory. Euphytica 47:189–194

Ikmal AM, Z Nurasyikin, A Kumar, AAS Noraziyah (2018). Evaluation of morpho-physiological traits of MRQ74 pyramided lines with drought yield QTLs. Euphytica 214:98

Ikmal AM, Z Nurasyikin, TATNA Riana, ZPD Ellina, R Wickneswari, AAS Noraziyah (2019). Drought yield QTL (QDTY) with consistent effects on morphological and agronomical traits of two populations of new rice (Oryza sativa) lines. Plants 8; Article 186

IRRI (2013). Standard Evaluation System for Rice. IRRI, Los Banos, The Philippines

Kang MS (1988). A rank-sum method for selecting high-yielding, stable corn genotypes. Cer Res Commun 16:113–115

Magari R, MS Kang (1993). “Genotype Selection via a New Yield-Stability Statistic in Maize Yield Trials.” Euphytica 70:105–111

Malaysian Meteorological Department (2017). Ringkasan Iklim 2016. Kuala Lumpur, Malaysia

Nassar R, M Huhn (1987). Studies on estimation of phenotypic stability: Tests of significance for nonparametric measures of phenotypic stability. Biometrics 43:45-53

Oladosu Y, MY Rafii, N Abdullah, U Magaji, G Miah, G Hussin, A Ramli (2017). Genotype × environment interaction and stability analyses of yield and yield components of established and mutant rice genotypes tested in multiple locations in Malaysia. Acta Agric Scand Sect B Soil Plant Sci 67:590–606

Othman S, ZPMD Hussain, CC Sheng, A Man, R Ahmad, B Amzah (2014). Padi aerob untuk mengatasi masalah kekurangan air. J Teknol 70:65–68

Poli Y, D Balakrishnan, D Subrahmanyam, M Panigrahy, SR Voleti, SK Mangrauthia, S Neelamraju (2018). Genotype × environment interactions of Nagina22 rice mutants for yield traits under low phosphorus, water limited and normal irrigated conditions. Sci Rep 8; Article 15530

Shahriari Z, B Heidari, A Dadkhodaie (2018). Dissection of genotype × environment interactions for mucilage and seed yield in Plantago species: Application of AMMI and GGE biplot analyses. PLoS One 13; Article e0196095

Shamsudin NAA, BPM Swamy, W Ratnam, MTS Cruz, A Raman, A Kumar (2016a). Marker assisted pyramiding of drought yield QTLs into a popular Malaysian rice cultivar, MR219. BMC Genet 17; Article 30

Shamsudin NAA, BPM Swamy, W Ratnam, MTS Cruz, N Sandhu, A Raman, A Kumar (2016b). Pyramiding of drought yield qtls into a high quality Malaysian rice cultivar MRQ74 improves yield under reproductive stage drought. Rice 9; Article 21

Shukla GK (1972). Some statistical aspects of partitioning genotype-environmental components of variability. Heredity 29:237–245

Sserumaga JP, SO Oikeh, S Mugo, G Asea, M Otim, Y Beyene, G Abalo, J Kikafunda (2015). Genotype by environment interactions and agronomic performance of doubled haploids testcross maize (Zea mays L.) hybrids. Euphytica 207:353–365

Suhaimee S, MR Rabu, AK Ali, A Ramli (2009). Ex ante assessment on feasibility of Maswangi (MRQ74) commercialization. Econ Technol Manage Rev 4:49–56

Swamy BPM, NAA Shamsudin, SNA Rahman, R Mauleon, W Ratnam, MTS Cruz, A Kumar (2017). Association mapping of yield and yield-related traits under reproductive stage drought stress in rice (Oryza sativa L.). Rice 10; Article 21

Thennarasu K (1995). On Certain Non-parametric Procedures for Studying Genoty Environment Interactions and Yield Stability. Indian Agricultural Research Institute, New Delhi, India

Tietenberg TH, L Lewis (2009). Environmental & Natural Resource Economics. Pearson, Boston, Massachusetts, USA

Vaezi B, A Pour-Aboughadareh, A Mehraban, T Hossein-Pour, R Mohammadi, M Armion, M Dorri (2018). The use of parametric and non-parametric measures for selecting stable and adapted barley lines. Arch Agron Soil Sci 64:597–611

Venuprasad R, CO Dalid, M Valle, D Zhao, M Espiritu, MTS Cruz (2009). Identification and characterization of large-effect quantitative trait loci for grain yield under lowland drought stress in rice using bulk-segregant analysis. Theor Appl Genet 120:177–190

Wricke G (1962). Uber eine methode zur erfassung der oekologischen streubreite in feldversuchen. Zeitschr Pflanzenz 47:92–96

Yan W, I Rajcan (2002). Biplot analysis of test sites and trait relations of soybean in Ontario. Crop Sci 42:11–20

Yan W, MS Kang (2002). GGE Biplot Analysis: A Graphical Tool for Breeders, Geneticists, and Agronomists. CRC Press, New York, USA

Yan W, MS Kang, B Ma, S Woods, PL Cornelius (2007). GGE biplot vs. AMMI analysis of genotype-by-environment data. Crop Sci 47:641–653

Yan W, LA Hunt, Q Sheng, Z Szlavnics (2000). Cultivar evaluation and mega-environment investigation based on the GGE biplot. Crop Sci 40:597–605

Zainuddin H, A Mokhtar, B Amzah, M Hashim, MBA Ghafar (2012). Enam varieti padi popular MARDI. Bull Teknol 1:1–10

Zhang H, Y Xue, Z Wang, J Yang, J Zhang (2009). Morphological and physiological traits of roots and their relationships with shoot growth in ‘super’ rice. Field Crops Res 113:31–40